function [error] = Modelfit_5d_parfor_toolbox(patient,ii,jj,kk,save_files);
    %| 5D model code - using the parallel toolbox
    
    %------------------------------------------------------------------------
    %|      Dependancies;                                                     |
    %|                  ParforProgMon (parfor loop progress monitor).
    %|                  read_dvf_elastix_slice
    %|                  iSaveX (parfor loop save function).
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 2014/04/01 10:23:59 $	$Revision: 0.1 $
    
    
    error=zeros(size(ii,2),size(jj,2),size(kk,2));%dim);
    % error_scans = zeros([size(run_scans,2), size(ii,2),size(jj,2),size(kk,2)]);
%     pctRunOnAll javaaddpath E:\ParforProgMonv2\java;% java path of ParforProgMon;
%     progressStepSize = 1;
%     ppm = ParforProgMon('          Lung: ',size(ii,2),  progressStepSize,1000, 100);
    
    tic
    
%     pctRunOnAll  warning('off')
    
    [~,grid_z]=meshgrid(jj,kk);
    grid_z=permute(repmat(grid_z,[1,1,patient.scans]),[3,2,1]);
    
    % block_size=8;
    % for block=1:round(size(ii,2)/block_size)
    
    for slice=ii;%((block-1)*block_size+1:(block-1)*block_size+block_size);
       
        if mod(slice, 10)==0
        display(sprintf('Fitting Slice %d of %d (%.0f %% done)', slice, size(ii,2), 100*slice/size(ii,2)))
        end
        [ dvf]=read_dvf_elastix_slice_toolbox(patient.dvf_folder,slice);
        
        dvf.u(patient.ref,:,:,:)=0;% Reference-to-reference DVF = 0;
        dvf.v(patient.ref,:,:,:)=0;
        dvf.w(patient.ref,:,:,:)=0;
        
        %%
        %Deform Bellows Voltage by dvf.w;
        X1=zeros(size(patient.run_scans,2),size(jj,2)*size(kk,2));%dim(2)*dim(3));
        X2=zeros(size(patient.run_scans,2),size(jj,2)*size(kk,2));%dim(2)*dim(3));
%         X3=zeros(size(run_scans,2),size(jj,2)*size(kk,2));%
                                                        %X3 = Cardiac term
        
        %Deform the bellows voltage index by the dvf.w;                                               
        voxel_w=double(bsxfun(@plus,grid_z,dvf.w(1:patient.scans,jj,kk)));
        voxel_w=voxel_w(1:patient.scans,:);
       
           for scan=1:patient.scans;
            X1(scan,:) = nakeinterp1([1:patient.dim(3)]',patient.bellows_volt_drifted(:,scan),voxel_w(scan,:)');
            X2(scan,:) = nakeinterp1([1:patient.dim(3)]',patient.flow_drifted(:,scan),voxel_w(scan,:)');
            %      X3(scan,:) = nakeinterp1(kk',ecg_volt(kk,scan),voxel_w(scan,:)');

        end
        %Vectorize the X1, X2 and DVF variables;
        X1=X1(patient.run_scans,:,:);X1=X1(:,:);
        X2=X2(patient.run_scans,:,:);X2=X2(:,:);
        
        dvf_u_vec=single(dvf.u(patient.run_scans,jj,kk));dvf_u_vec=dvf_u_vec(:,:);
        dvf_v_vec=single(dvf.v(patient.run_scans,jj,kk));dvf_v_vec=dvf_v_vec(:,:);
        dvf_w_vec=single(dvf.w(patient.run_scans,jj,kk));dvf_w_vec=dvf_w_vec(:,:);
        
        %Initialize the model params;
        tempfit_variables_x=zeros(3,length(X1));
        tempfit_variables_y=zeros(3,length(X1));
        tempfit_variables_z=zeros(3,length(X1));
        
%         temp_array_x=( zeros(3,1,length(X1)));
%         temp_array_y=(zeros(3,1,length(X1)));
%         temp_array_z=(zeros(3,1,length(X1)));
        
        
        yfit_x=zeros( size(patient.run_scans,2),length(X1));
        yfit_y=zeros( size(patient.run_scans,2),length(X1));
        yfit_z=zeros( size(patient.run_scans,2),length(X1));
        
        constant=zeros(length(jj),length(kk),3);
        alpha=zeros(length(jj),length(kk),3);
        beta=zeros(length(jj),length(kk),3);
        % gamma=zeros(length(jj),length(kk),3);
        
        
        %% For loop model fitting code;
%         tic
        for vox=1:length(X1);
            
            E_vox=[ones(size(X1,1),1) X1(:,vox) X2(:,vox)];
            
            tempfit_variables_x(:,vox)=E_vox\dvf_u_vec(:,vox);%(E_vox'*E_vox)\(E_vox'*dvf_u_vec(:,vox));%
            tempfit_variables_y(:,vox)=E_vox\dvf_v_vec(:,vox);%(E_vox'*E_vox)\(E_vox'*dvf_v_vec(:,vox));%
            tempfit_variables_z(:,vox)=E_vox\dvf_w_vec(:,vox);%(E_vox'*E_vox)\(E_vox'*dvf_w_vec(:,vox));%
            
            yfit_x(:,vox)=tempfit_variables_x(:,vox)'*E_vox';
            yfit_y(:,vox)=tempfit_variables_y(:,vox)'*E_vox';
            yfit_z(:,vox)=tempfit_variables_z(:,vox)'*E_vox';
            
        end
%         time_for=toc
                
        error_u = bsxfun(@minus,yfit_x,dvf_u_vec);
        error_v = bsxfun(@minus,yfit_y,dvf_v_vec);
        error_w = bsxfun(@minus,yfit_z,dvf_w_vec);
        error_tmp_scans=(bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v)));
        error_tmp=mean(bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v)));
        error(slice,:,:)=reshape(error_tmp,size(jj,2),size(kk,2));
        
        
        if save_files==1
                        
            constant(jj,kk,1)=reshape(tempfit_variables_x(1,:),size(jj,2),size(kk,2));
            constant(jj,kk,2)=reshape(tempfit_variables_y(1,:),size(jj,2),size(kk,2));
            constant(jj,kk,3)=reshape(tempfit_variables_z(1,:),size(jj,2),size(kk,2));
            
            alpha(jj,kk,1)=reshape(tempfit_variables_x(2,:),size(jj,2),size(kk,2));
            alpha(jj,kk,2)=reshape(tempfit_variables_y(2,:),size(jj,2),size(kk,2));
            alpha(jj,kk,3)=reshape(tempfit_variables_z(2,:),size(jj,2),size(kk,2));
            
            beta(jj,kk,1)=reshape(tempfit_variables_x(3,:),size(jj,2),size(kk,2));
            beta(jj,kk,2)=reshape(tempfit_variables_y(3,:),size(jj,2),size(kk,2));
            beta(jj,kk,3)=reshape(tempfit_variables_z(3,:),size(jj,2),size(kk,2));
            
            iSaveX( [patient.model_params_folder sprintf( '/constant_%d', slice )], constant );
            iSaveX( [patient.model_params_folder sprintf( '/alpha_%d', slice )], alpha );
            iSaveX( [patient.model_params_folder sprintf( '/beta_%d', slice )], beta );
        end
        % toc
        
%         ppm.increment();
        
    end
    
    if save_files==1
        save([patient.model_folder '/error'], 'error')
        
    end
    
    alpha = compile_3d_alpha_toolbox(patient)
    beta = compile_3d_beta_toolbox(patient)
    constant = compile_3d_constant_toolbox(patient)

    toc
    
    